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Summary. We present a method to decompose an arbitrary 3D piecewise linear 
complex (PLC) into a constrained Delaunay tetrahedralization (CDT). It success- 
fully resolves the problem of non-existence of a CDT by updating the input PLC into 
another PLC which is topologically and geometrically equivalent to the original one 
and does have a CDT. Based on a strong CDT existence condition, the redefinition is 
done by a segment splitting and vertex perturbation. Once the CDT exists, a prac- 
tically fast cavity retetrahedralization algorithm recovers the missing facets. This 
method has been implemented and tested through various examples. In practice, it 
behaves rather robust and efficient for relatively complicated 3D domains. 



1 Introduction 

A fundamental problem in unstructured mesh generation is to create a mesh that 
represents a geometric domain bounded by piecewise linear faces and possibly with 
interior constraining faces. It is also referred as boundary mesh generation. Numer- 
ous applications depend on it. 

This problem has been successfully solved in two dimensions. It is well known 
that every polygonal domain can be triangulated into triangles without adding new 
vertices (the Steiner points). Almost optimal algorithms (with a linear complexity 
in practice) have been proposed [Lee86, Chew89]. 

The problem is significantly more difficult for arbitrarily shaped three dimen- 
sional domains. Such domains can be described by piecewise linear complexes 
(PLCs) [Miller96] , which are objects more general than polyhedra. It is known [Schoen- 
hardt28, Bagemihl48, Chazelle84, Rambau03] even a simple polyhedron may not be 
tetrahedralizable without adding new vertices. The simplest example is the so-called 
Schonhardt polyhedron [Schoenhardt28] , which is a non-convex twisted triangular 
prism. Moreover, the problem of deciding whether a simple polyhedron can be tetra- 
hedralized is NP-hard [Ruppert92]. PLCs are usually much more complicated than 
simple polyhedra. To guarantee an arbitrary PLC can always be meshed, methods 
must resort to adding Steiner points. However, a number of difficult issues remain 
to be resolved, such as the placement of the Steiner points, the minimum bound on 
such points, and so on. 
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The recently proposed Constrained Delaunay tetrahedralization (CDT) (cl. [Shew- 
chuk02]) is a Delaunay-like tetrahedralization that is constrained to respect the 
shape of a PLC. CDTs are obviously suitable structures for resolving the above 
problem. Not only they respect the boundary but also they have many nice mathe- 
matical properties inherited from Delaunay tetrahedralizations. Many applications 
can be envisaged after getting a CDT. For instance, it is a good initial mesh for 
getting a quality conforming Delaunay mesh [Shewchuk98b, Cheng04, Pav04] which is 
suitable for numerical methods. 

A key question for constructing a CDT is to decide its existence, i.e., whether 
a given PLC has a CDT without adding points. So far, Shewchuk [Shewchuk98a] 
has proved the condition: if all segments of the PLC are strongly Delaunay, then 
the CDT exists. The hint gained from this condition is: additional points can be 
inserted only on the segments of the PLC. 

Shewchuk [Shewchuk02] gave a segment recovery algorithm for constructing 
CDTs. For a PLC X, it carefully introduces additional points on some segments of 
X until the existence of a CDT can be guaranteed. This algorithm yields a provably 
good bound on edge lengths. However, it requires to compute the local feature size 
explicitly for protecting sharp corners (vertices with angles less than 90° formed by 
segments). In [Si04a] we proposed a new segment recovery strategy which exploits 
the available local geometric information to efficiently construct Steiner points on 
segments. It needs not to compute the local feature sizes. Moreover sharp corners 
are implicitly handled during the creation of the Steiner points. Both algorithms 
tend to use fewer additional points than other methods which do edge protect prov- 
ably [Pebay98, MurphyOO, Cohen-Steiner02] too. 

When the existence of a CDT is known, another key issue is to recover facets 
of the PLC. Generally non-CDT algorithms [Pebay98, MurphyOO, Cohen-Steiner02, 
Weatherill94, George03, Du04] will continuously insert points on the missing facets 
or the inside of the PLC. While CDT algorithms [Shewchuk02, Shewchuk03, Si04a] 
recover the missing facets without introducing additional points. This again reduces 
the number of Steiner points in a CDT. 

By now Shewchuk has provided several facet recovery algorithms [ShewchukOO, 
Shewchuk02, Shewchuk03]. The incremental facet insertion algorithm [Shewchuk02] 
recovers facets one after one and the CDT is updated accordingly. For each facet, a 
gift- wrapping algorithm is used for retetrahedralizing the two cavities around it. The 
strategy is simple but the time complexity is poor and the performance is unstable 
due to the gift-wrapping algorithm. The flip-based algorithm [Shewchuk03] recovers 
each facet by a sequence of carefully ordered flip operations. It appears to be simple 
and is likely to outperform other algorithms on most inputs. In order to guarantee 
the correctness and termination, both algorithms require that a full perturbation 
has to be applied on the set of vertices to remove the degeneracies. 

In this paper we present a new CDT algorithm. The main difference from other 
CDT algorithms is the practical exploitability of a strong CDT existence condition 
which requires no local degeneracy on the vertices of the PLC. We propose a local 
degeneracy removal algorithm to construct a new set of vertices out of the old one 
which is consistent with the constraining segments and facets of the PLC. After the 
strong condition is satisfied, facet recovery is done by a new cavity retetrahedraliza- 
tion algorithm which is fast and robust in practice. 

The remainder of this paper is organized as follows. We shortly recall the defini- 
tions of PLCs and CDTs in the next section. Then the existence of CDT is discussed 
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in section 3. Section 4 provides an overview of the proposed method. The individ- 
ual algorithms are fully described in section 5, 6, and 7. Finally we present some 
meshing results from publicly available examples. 



2 PLCs and CDTs 

A piecewise linear complexes (PLC) proposed by Miller, Teng, Walkington, and 
Wang [Miller96] is a general boundary description for three-dimensional domains. 
Simply saying, a PLC X is a set of vertices, together with a collection of segments, 
and facets. Like a simplical complex, any two components of X are either disjoint or 
meet in a common face, e.g., two segments can only intersect at a vertex of X, and 
two facets can only intersect at a collection of segments and vertices of X, and so 
on. A PLC facet has no analogue in a simplical complex. Each facet of X is indeed 
a two-dimensional polygonal region embedded in three dimensions, it may not be 
convex and possibly contains holes, isolated segments and vertices. 

PLCs are more general than polyhedra in the sense that every polyhedron is a 
PLC but not vice versa. For instance, a PLC containing a segment inside can't be 
represented by any polyhedron. Surface triangulations are one special type of PLCs - 
each facet is a triangle. Hence PLCs are able to approximate arbitrary complicated 
and curved shapes. In addition, many popular polygonal file formats (e.g., STL, 
OFF, PLY, and [Si04b]) can be directly used or slightly modified to describe PLCs. 

Given a PLC X, Shewchuk [Shewchuk02] defined a CDT of X as follows: 

Let V be the set of vertices of X. a is any simplex (tetrahedron, triangle, edge 
or vertex) formed by vertices of V . a is Delaunay if there exists a circumsphere 
of a that encloses no vertex of V . The Delaunay tetrahedralization T> of V is a 
tetrahedralization that all simplices of T> are Delaunay. T> is unique if V is general, 
i.e., no five or more vertices lie on a common sphere. 

The visibility between two vertices p and q is occluded if there is a constraining 
facet / such that p and q lie on opposite sides of the plane that includes /, and the 
line segment pq intersects this facet (see Figure 1). Segments do not occlude visibility. 
For example, in Figure 1, c and d can see each other even if ab is a segment. 

A tetrahedron t formed by vertices of X is constrained Delaunay if its circum- 
sphere encloses no vertex of X which is visible from any point in the relative interior 
of t (see Figure 1). 

A tetrahedralization T is a constrained tetrahedralization of X if T and X have 
exactly the same vertices, and the tetrahedra in T cover the convex hull of V. Every 
segment of X is an edge of T, every facet of A is a union of triangular faces of T . 

A constrained tetrahedralization T of X is said to be a constrained Delaunay 
tetrahedralization (CDT) of X if each tetrahedron of T is constrained Delaunay. 

Intuitively, the definitions of Delaunay tetrahedralization and constrained De- 
launay tetrahedralization are the same except that, for the CDT, we ignore the 
volume of a sphere whenever the sphere passes through a facet of X. 

Let T be a CDT of A. A facet of X is represented by a set of coplanar triangular 
faces of T. Such faces are called subfaces for distinguishing them from other faces 
of T. The set of subfaces of a facet form a two-dimensional constrained Delaunay 
triangulation. 
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Fig. 1. Constrained Delaunay tetrahedron. The shaded region represents a facet / 
of X including vertices a, b, c and d. Vertices p and q lie on opposite sides of /, they 
can not see each other. While c and d can see each other even if ab is a segment of 
X. St is the circumsphere of tetrahedron t (abcp) and it encloses q but not d, t is 
constrained Delaunay. 

3 The Existence of CDT 

Given a PLC X. Generally, the CDT of X may not exist. One reason is that X may 
not be tetrahedralizable at all. Even the constrained tetrahedralization of X exists, 
it may still not have a CDT. A key question for CDT algorithms is to decide under 
what condition the CDT exists. 

Shewchuk has proved a condition. Let a be a simplex of X, a is strongly Delaunay 
if there exists a circumsphere of a that passes through and encloses no other vertices 
of X. Say X is edge-protected if every segment of X is strongly Delaunay. 

Theorem 1 ( [Shewchuk98a]). If X is edge-protected, then X has a CDT. ■ 

This condition is useful in practice because it suggests that additional points can 
be inserted on segments only. However, this condition does not help us to construct 
the CDT. Suppose X has been made edge-protected (by splitting the segments into 
smaller segments), algorithms [Shewchuk02, Shewchuk03, Si04a] need to carefully 
do perturbations on the vertices of X in order to successfully recover the missing 
subfaces. 

In this paper, we use another less general condition which guarantees the exis- 
tence of CDT, too. 

Definition Let T> be the Delaunay tetrahedralization of the vertices of X, t a 
tetrahedron in T> and t' an adjacent tetrahedron of t (sharing a face with t), V the 
set of vertices of t, t' . If all vertices of V lie on a common sphere, V is called a local 
degeneracy of X. 

Remark. If X contains no local degeneracy, then T> is unique. However, the set 
of vertices of X may still be degenerate ( since arbitrary 5 or more vertices of X 
can lie on a common sphere). 

Theorem 2. IfT> contains no local degeneracy and contains all segments of X , then 
the CDT of X exists. 
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Proof: It is easy to show: if X satisfies the above condition, then the segments of 
X are edge-protected. ■ 

This condition is stronger than Shewchuk's condition since it implicitly satisfies 
it. Its advantage is: the facet recovery can be carried out efficiently. Section 7 presents 
such an algorithm. 



4 The Algorithm 

Let the initial PLC be Xo. The CDT is constructed by the following consecutive 
steps: 

(1) Construct an initial Delaunay tetrahedralization 2?o of the vertices of Xo. 

(2) Recover the segments of Xo in T>o by incrementally inserting points on missing 
segments, update Xo —> X\ and Do — ► T>\ with the newly inserted points 
respectively. 

(3) Remove the local degeneracies in X\ by either perturbing vertices or inserting 
new vertices, update X\ — > X2 and 27 1 — > 2?2 with the newly inserted points 
respectively. 

(4) Recover the subfaces of X2 in T>2 by a cavity retetrahedralization method. 

In step (1), 2?o can be efficiently constructed by any standard algorithm, such 
as [Edelsbrunner96] . T>o probably does not respect the segments and subfaces of Xo- 
After step (2) is done, D\ is a Delaunay tetrahedralization of vertices of X\ which 
contains all segments of Xi . However, Di may not respect the subfaces of Xi . Step 
(3) guarantees the existence of a CDT. After all local degeneracies are removed, 2?2 
is the unique Delaunay tetrahedralization of vertices of X2 and contains all segments 
of X2. D2 may still not respect the subfaces of X2, while the existence of the CDT (of 
X2) is guaranteed. Hence the step (4) can be carried out without inserting additional 
vertices. These processes are detailed in the following sections. 



5 Segment Recovery 

A segment of Xo is missing if it is not in 2?o- The purpose of the segment recovery 
algorithm is to update T>o into 2?i such that 2?i is a Delaunay tetrahedralization 
and includes all segments of Xo- 

Let dej be a segment with endpoints e; and ej, be its Euclidean length. 

A vertex is acute if at least two segments incident at it form an angle smaller than 
90°. We distinguish two types of segments, a segment is type-1 if its both endpoints 
are not acute, it is type-2 if only one of its endpoints is acute. If both endpoints of 
a segment are acute, it can be transformed into two type-2 segments by inserting a 
vertex. 

A segment is split into subsegments. Subsegments inherit types from original 
segments. For example, let aej be a subsegment of eie2 which is a type-2 segment 
and ei is acute, is type-2 although none of its endpoints is acute. For any vertex 
v inserted on a type-2 segment (or subsegment), O(v) denotes its original acute 
vertex. A tacit rule is used throughout this section, if eiej is type-2, it implies either 
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ei or O(ei) is the acute vertex. In the following, unless it is explicitly mentioned, a 
segment can be a segment or subsegment. 

A vertex encroaches upon a segment if it lies inside or on the diameter sphere 
of that segment. Remark: any missing segment must be encroached by at least one 
vertex of X. 

Let dej be a missing segment, it will be split by a vertex v. The reference point 
p of v, which is "responsible" for the insertion of v, is defined as follows: 

• p encroaches upon e^ ; 

• the circumradius of the smallest circumsphere of triangle dejp is maximum over 
other encroaching points of e;e-j. 

Figure 2 illustrates how p is chosen. Notice that p may not unique (because several 
points can share the same circumsphere), choose an arbitrary one in this case. 



-Pa 

P2 • ' 



• Pa = P . 
e< *f % e. 



Fig. 2. The reference point p of a splitting segment etej. pi, p2, P3 and p4 all 
encroach upon etej. p4 is chosen as the reference point because it forms the biggest 
circumsphere with e^ej. 



The insertion of v to split eiej is governed by three rules given below. Let S(c, r) 
denote a sphere centered at c with radius r: 

1. dej is type-1 (see Figure 3 (a)), then v := eiej n S, where S(c,r) is the sphere 
defined by the reference point p of v as follows: 

if \eip\ < i|e»ej| then 

c := e t ,r := |e;p], 
else if \ejp\ < ||e,:ej| then 

c := ej,r:= \e jP \, 
else 

c . — ei,r . — 2!^*^ I' 
end 

2. eiej is type-2 (see Figure 3 (b)), let eu := O(ei), then jj := e^ej n S, where 
S(c, r) is the sphere defined by the reference point p of v with c := et,r := |ejtp| . 
However, if the subsegment vej has length \vej\ < \vp\, then we do not insert v 
and use rule 3 to split the segment. 
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3. eiej is type-2, let e k := 0(ei), and v' is the rejected vertex by rule 2. then 
v := e^ej n S, where S(c,r) is the sphere defined by the reference point p of v 
as follows (see Figure 3 (c)): 
c := e k 

if \pv'\ < §|eii/| then 

r := |efcfii | + \en>'\ - \pv'\ 
else 



end 



r := \e k ei\ + ||ejv'| 



(a) (b) (c) 



Fig. 3. Illustrations of the segment splitting rules. 

For several segments sharing an acute vertex, by repeatedly using rule 2 or 3, 
a protecting ball is automatically created which ensures: no other vertex can be 
inserted inside the ball. The effect is shown in Figure 4. Notice, the protecting ball 
is not necessarily completely created, only the missing segments will be split and 
protected. Existing segments remain untouched. This reduces the number of Steiner 
points. 

Below is the pseudo-code of the segment recovery algorithm. 

Algorithm Delaunay Segments Recovery 
Input: T>o, Xq. 
Output: 2?i, Xl 
initialize: 

T>\ := T>o, Xi := Xo; 
repeat: 

form a queue Q of missing segments in T)\ ; 
while Q j= 0 do 

remove a segment e^ej from Q; 

split eiej using rule 1, or 2, or 3; 

update T>i, Xi; 
end 

until no segment of X\ is missing in T>\ 
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Fig. 4. The protecting ball of an acute vertex e*. Vi, V2, i>3, and, v± are points 
inserted on segments (by rule 2) sharing e». They automatically create a protecting 
ball of e t . 

The termination of this algorithm can be proved by showing that the length 
of every segment created is bounded by the local feature size divided by constant 
depending only on the input. For a PLC X, the local feature size [Ruppert95] lfs(v) 
of any point v in X is the radius of the smallest ball centered at v that intersects 
two segments or vertices in X that do not intersect each other. The lfs() defines 
a continuous map that maps every point in X into a positive value which suggests 
how large the ball of the empty space around this point can be. The function lfs() 
is only defined on the input PLC X and does not change as new points are inserted. 

Theorem 3 ([Si04a]). Let eiej be a finally resulting subsegment, 

• if etej is type-1, then: 

|e»ej| > min{lfs(ei),lfs(ej)}. 

• if eiej is type- 2, let e k := 0(d), then: 
|e»ej| > ^lfs(ek) when e< = e k , 
\ e i e j \ > lfs(ek)sin(6) when a / e k . 

where C is two times the number of segments incident at e k and 9 is the smallest 
angle between them. 

hence the Delaunay segments recovery algorithm terminates. ■ 

Practically, the algorithm terminates within a few steps creating the protection 
ball sector. The constant C usually is no larger than 4. 



6 Removing Local Degeneracies 

In order to fulfill the condition of Theorem 2, all local degeneracies have to be re- 
moved from 2?i . Techniques of perturbation [Edelsbrunner90] are effective to remove 
degeneracies. However, they must be carefully applied in CDT algorithms. 

We say a vertex of a PLC is perturbable if there exists an arbitrarily small 
perturbation on it which does not affect the consistency of the PLC, otherwise, it is 
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imperturbable. For example, a vertex on a segment or inside a facet is perturbable 
since it can be perturbed arbitrarily along the segment vector or inside the facet 
without affecting the collinearity of the segment or the coplanarity of the facet. Not 
every vertex of a PLC is simply perturbable. If a vertex intersected by three or more 
non-coplanar facets is perturbed, at least one facet becomes invalid (because it now 
contains a non-coplanar vertex), hence it is imperturbable. 

Let A be a local degenerate set of vertices in Xi, i.e., A contains 5 vertices of 
Xi which share a common sphere. If p £ A is perturbable, A can be removed by an 
arbitrarily small perturbation on p. We call a perturbation is segment-safe if after 
the perturbation no segment of X\ becomes non-Delaunay. A perturbable vertex 
may not be segment-safe, see Figure 5 for an example. 




Fig. 5. A perturbable vertex which is not segment-safe. The vertices of this facet 
consist of a 3 x 3 square grid. The vertex p is perturbable but not segment-safe when 
the four edges opposite it are all segments. Whatever it moves in the facet, at least 
one segment becomes non-Delaunay. 

A is said to be removable if there exits a vertex p € A, such that: 

• p is perturbable; and 

• there exists a perturbation of p which is segment-safe. 

Otherwise, A is unremovable. 

If A is unremovable. We introduce a new point vt, called break point, it is chosen 
as follows: 

• If the vertices of A are affmely independent. Let Sa be the common sphere 
shared by them, vt, is inside Sa- 

• If the vertices of A are affinly dependent, i.e., four of them are coplanar. Let 
Ca be the common circle shared by the four vertices, vt is coplanar with the 
four vertices and inside Ca- 

Vb will break the local degeneracy (see Figure 6). 

A break point may encroach upon one or more segments and subfaces of Xi . In 
this case, it will not be inserted. Instead, the following boundary protection procedure 
is called: 

1. For each encroached segment, add its perturbed circumcenter to Xi- Updating 
2?2, X2 accordingly; 
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Fig. 6. The break point. On the left is a set of locally degenerate points and one of 
its Delaunay triangulations. Vb is a break point for removing the local degeneracy. 
On the right is the unique Delaunay triangulation after vt is inserted. 

2. For each encroached subface, compute its perturbed circumcenter x. If x en- 
croaches upon any segments, it is not inserted, goto step 1 to split the encroached 
segments. Otherwise, add x to X2. Updating T>2, X2 accordingly; 

3 Call the Delaunay segment recovery algorithm to recover all missing segments 
of X2 in V2 ■ 

The complete algorithm is listed below: 

Algorithm Local Degeneracy Removal 
Input: 2?i, X±. 
Output: V 2 , X 2 . 
initialize: 

V 2 :=T>i, X 2 ~Xu 
repeat: 

form a queue Q of local degeneracies of T>2 ; 
while Q / 0 do 

remove a local degeneracy A from Q; 

if A is removable then 

Remove A by a small perturbation; 

else 

Compute a vi of A; 

if Vb encroaches upon any segment and subface then 
Push A into Q; 

Call the boundary protection procedure; 
else 

Insert Vb to break A; 
update T>2, X2; 
endif 
endif 
endwhile 
until T>2 contains no local degeneracy; 

In our implementation, we choose Vb be the circumcenter of the common sphere 
of A. It is important that Vb should not create a new local degeneracy with other 
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existing vertices. This can be checked before inserting the point. If so, we should 
not insert vi but choose another location. Notice that such case only can happen 
when the input PLC is highly symmetric, e.g., in Figure 5. A simple strategy is to 
add a randomized perturbation on each «(,. After the perturbation, the probability 
to create a symmetric point configuration again is nearly zero. Due to that reason, 
the points added for boundary protecting are perturbed, too. 

Some break points may locate outside X2. They will be removed after the CDT 
of X2 is constructed. 

The termination of the local degeneracy removal algorithm is due to the fact 
that the number of unremovable local degeneracies can only be decreased and no 
new local degeneracies are created by break points and segment protecting points. 



7 Facet Recovery 

The segment recovery and local degeneracy removal algorithms have produced X2 
such that it has a CDT (guaranteed by Theorem 2). Let T be a CDT of X2. Gen- 
erally, T>2 7^ T because some subfaces of T are non-Delaunay faces and penetrated 
by edges of T>2- This section describes an algorithm which incrementally transforms 
T>2 into T. No additional points are needed in the transformation. 

The general steps of our algorithm are similar to [Shewchuk02] . At initialization, 
let T (0) := £> 2 ; add all missing subfaces into a queue Q. The algorithm starts to 
recover the subfaces in Q until Q is empty. At each step i, the algorithm recovers a 
set of missing subfaces in update into T (i+1) . After m steps Q is empty, 

and T = T (m) . 

At step i(i < m), several missing subfaces are recovered together. We define a 
missing region J? to be a set of coplanar subfaces of X2 such that 

• all subfaces of Q belong to one facet of X2; 

• the boundary edges of Q are edges of T' 1 '; and 

• the internal edges of Q are missing in T*- 1 ' . 

Hence Q is a connected set of missing coplanar subfaces. It may not simply con- 
nected, i.e., Q can contain a hole inside (see Figure 7 (a)). Each missing subface 
belongs to one missing region. A facet can have more than one missing regions. 



Fig. 7. (a) The shaded area highlights a non-simply connected missing region Q. 
(b) A cavity C at (step i) is illustrated. 




(a) 



(b) 
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When a Q is found, the formcavity procedure (described below) forms two cav- 
ities at each side of Q. Each cavity is bounded by subfaces of fl and faces of 
(see Figure 7 (b)). 

1. Find all tetrahedra in that intersect the relative interior of Q, delete them 
from T^' . This creates a hole inside . 

2. Insert the missing subfaces of Q into the hole to split it into two separated 
cavities, one at each side of O. 

Each cavity C is filled with a set of new tetrahedra by the following cavity 
retetrahedralization procedure. Let V be the set of vertices of C . 

1. Verify C, expand C if it is necessary. 

(1) Form a queue Q containing all non-strongly Delaunay faces of C in V; 

(2) For each face a £ Q and a still in C, 

let t be the tetrahedron adjacent to C and holds a; 

remove a from C; 

for each face At of t, At ^ a do 

if At is a face of C then remove At from C; 

else add At into C; end 
end 

Update C, V; 

(3) Repeat (1) if some faces of C are not strongly Delaunay in V . 
See Figure 8 for an example of the expansion of C. 

2. Retetrahedralize C. 

(1) Construct a Delaunay tetrahedralization T>c of the vertices of the C. 

(2) Identify faces of C in T>c, mark each tetrahedron of T>c to be "inside" or 
"outside" . 

(3) Remove tetrahedra marked as "outside" from T>c, and fill the remaining 
tetrahedra into C. 

Figure 9 illustrates a two dimensional example of the retetrahedralization pro- 
cedure. 




Fig. 8. The expansion of the cavity (illustrated in 2D). Left: eie2 is the segment 
going to recover. Edges below eie2 is the faces of the cavity C. piq2 is not strongly 
Delaunay. Right: C is expanded by removing two edges piqi and piq2 from C and 
adding one edge 51^2 into C. 
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(a) (b) (c) (d) 

Fig. 9. Cavity triangulation (illustrated in 2D), (a) Two initial cavities separated 
by a constraining segment, (b) The two Delaunay triangulations constructed at 
each side of the segment, (c) Mark triangles as "inside" or "outside", (d) Remove 
"outside" triangles. 

The complete facet recovery algorithm is summaried as follows: 

Algorithm Facet Recovery 
Input: 2?2, X%. 
Output: The CDT T of X 2 . 
initialize: 

T:=T> 2 ; 
Repeat: 

form a queue Q of missing subfaces in T; 
while Q / 0 do 

remove an unrecovered subface / from Q; 

form a missing region Q containing /; 

form two cavities Ci, C2 by formcavity subroutine; 

for each d, i = {1, 2} do 

Call cavity retetrahedralization subroutine; 

end 
end 

until no subfaces are missing in T. 
Theorem 4. The facet recovery algorithm terminates. 

Proof: We show that the cavity retetrahedralization subroutine will always succeed, 
hence the missing region found at each step can be recovered without getting stuck. 

The step 1 (face verification) of the cavity retetrahedralization algorithm guar- 
antees the step 2 can be correctly executed since every strongly Delaunay simplex 
of V will appear in any Delaunay tetrahedralization of V . Hence the face identifi- 
cation in step 2 must be successful. What remains is to prove two issues in the face 
verification step: (1) the expansion of C terminates; and (2) the missing subfaces 
due to the expansion of C can be recovered later. 

Let 9 be the set of faces of T>i such that no face of ^ is crossed by any subfaces 
of X 2 . Clearly, •f' 7^ 0 and any face At € •f' is strongly Delaunay and exists in any 
T^'. The set <P limits the expansion of C, i.e., C stops expanding at a when 0 £ ^ ■ 

Let Vol(C) be the inside volume of C . During the expansion of C some subfaces 
are missing. Notice that the missing region Q' formed by these missing subfaces is 
completely inside C, hence Vol(C') < Vol(C), where C" is the cavity formed from 
Q' . This relation holds if the expansion of C' still causes some other subfaces miss- 
ing and results another new cavity. Such sequence will terminate since the value of 
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VolQ can not be negative. 



Theorem 5. T created by the facet recovery algorithm is a CDT. 

Proof: We first show T' 1 ' is a CDT. is the result of cavity retetrahedralization 
algorithm on T^°\ Tetrahedra of T^ 1 ' which are outside and not adjacent to the 
cavities remain Delaunay. Let t be a tetrahedron created inside a cavity C, St is its 
circumsphere. Let t' be another tetrahedron in T' 1 ' sharing a face a with t, v be 
the vertex of t! opposite to a. We have the following cases: 

(1) a is a subface. Then t is constrained Delaunay even if St encloses v, i.e, the 
inside of t is not visible by v, otherwise at least a segment is non-Delaunay and 
can not exist in T' 0 ' ; 

(2) a is a face inside C, then St must not enclose v (guaranteed by the cavity 
retetrahedralization algorithm) . 

(3) a is a face on C, then a is not a subface, St must not enclose v. Otherwise, T^ 0 ' 
is not a Delaunay tetrahedralization since the circumsphere of t' is not empty, 
i.e., it encloses the point of t opposite to a. 

By induction, after step i, i > 1, T' 1 ' is a CDT. Now we can show 7"'* +1 ' is a 
CDT by using the similar arguments as above. The only difference is in the case (3) 
which is stated below: 

(3) a is a face on C, then we have two cases, 

a. a is not a subface, St must not enclose v. Otherwise, T'*' is not a CDT 
since the circumsphere of t' encloses the point of t opposite to a which is 
also visible from the inside of t! . 

b. a is a subface, then t is constrained Delaunay even if St encloses v. 

Thus on the finish of the facet recovery algorithm T = T' m ' is a CDT. ■ 



8 Experimental Results 

This algorithm has been implemented and in our 3D quality Delaunay mesh gen- 
erator - TetGen [Si04b] (http://tetgen.berlios.de). We have tested the program 
with a number of examples not only having simple but also relatively complicated 
geometries. The algorithm runs rather efficiently compares to the old version of Tet- 
Gen which uses gift-wrapping algorithm. For example, for a cavity having around 
300 faces and 150 vertices, the gift-wrapping algorithm needs few minites to finish 
while the cavity retetrahedralization algorithm finishes in less than 1 second. To 
test the stability of the algorithm as well as our implementation, we purposely chose 
some models which the surface meshes are rather badly discretized. They are most 
likely to pull down the program if the algorithm has defect. On most of these models 
the program successfully produced CDTs as long as the surface meshes are PLCs. 

The following two PLC models can be freely downloaded from Inria's large repos- 
itory http: //www-rocql . inria. fr /gamma. 

The Cup-holder (in Figure 10) has a simple geometry (918 vertices, 1848 trian- 
gles). 1261 points are added (262 break points and 999 protecting points). There are 
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157 missing subfaces, the biggest size of a cavity contains 40 faces. Moreover, the 
effect of removing local degeneracies can be visualized on the right picture. 

The monster4 (in Figure 11) consists of 1392 vertices, 2784 triangles. It has 
complicated distribution of segments and facets as shown in the left picture. Our 
algorithm added 2782 Steiner points (502 break points and 2226 protecting points). 




Fig. 10. Example 1 (Cup-holder). Left: the PLC. Right: the CDT. 




Fig. 11. Example 2 (monster4). Left: the PLC. Right: the CDT. 

The geometry of Figure 12 is a human heart inside a body. The surface mesh 
of the heart is internal boundary separating two regions. Hence this model is made 
up of multiple domains. In spite of the complexities in the geometries, the point 
set contain no local degeneracy. Hence the CDT only require few additional points 
(compare to the input number of points) for protecting segments. The CDT is shown 
in the middle. One can see that the tetrahedra of the CDT are usually very long 
and skinny. The picture on the right shows the better quality tetrahedra obtained 
by performing a Delaunay refinement on the CDT. The number of additional points 
for improving the mesh quality are much bigger than that of getting a CDT. 
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Fig. 12. Example 3 (Heart). Left: the PLC. Middle: the CDT. Right: the quality 
mesh. A vertical cut has been made on the meshes so that the internal boundaries 
can be seen. 
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